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Splitting a solution of the wave equation in n + 1 space dimensions into spherical harmonic com- 
ponents on n-spheres and reducing to first order results in the pair of first-order partial differential 
equations in radius and time -k = ip'+pip/r and ip = vr' for each spherical harmonic, where p — 2l + n, 
and I is a spherical harmonic index. The ip/r term gives rise to numerical stability problems near 
the origin, and also poses the key numerical difficulty in related systems of equations. We propose a 
class of summation by parts finite differencing methods that converge pointwise, including at r = 0, 
and are stable because they admit a discrete energy. We explicitly construct such schemes that are 
2nd and 4th order accurate at interior points, and first and 2nd order accurate at the spherical outer 
boundary r — R. We use the projection method to impose a class of boundary conditions for which 
the discrete energy is non-increasing. 
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I. INTRODUCTION 

We consider the wave equation 

i> = a$ 



(1) 



in n + 1 space dimensions and time. In spherical coordi- 
nates, the Laplace operator A in n + 1 space dimensions 
can be split into radial and angular derivatives as 



A = 



J_d_ 

r n dr 



dr 



1 



A 



(2) 



The Laplacian on S n has eigenfunctions Yj... that obey 
A sn Y L .. = -l(l + n-l)Y L ... (3) 
We can therefore make the separation of variables ansatz 



$(r,M)=^J> ±i (r,i)e d 
z=o ± 

in 2 space dimensions, 

oo Z 



(4) 



(5) 



Z=0 m=-l 



in 3 space dimensions, and 



$(r,*, angles) = ^^<^...(r,i)YL. (angles) (6) 
i=o ... 

in higher space dimensions, where the Y; m are the famil- 
iar scalar spherical harmonics, and Yj... their generalisa- 
tion to higher dimension. The nonstandard notation for 
the Fourier series in the n = 1 case was chosen to give I 
the range I > for all n > 1 . n = 1 is then a regular case 
in what follows. Each </>/... now obeys 



1(1 + n- 1) 



(7) 



Here a dot denotes a ^-derivative, and a prime an r- 
derivative. We no longer write the suffix ±Z, ^to, or 
I . . . that labels the spherical harmonic. Regularity of 
$ at r = 0, in the sense that it admits an expansion 
in positive integer powers of the Cartesian coordinates 
x,y,z, requires that 4>(r,t) = r l cf>(r,t), where (j> is even 
and regular, in the sense that it admits an expansion in 



positive even powers of r. 



J ; 



is difficult to enforce 



numerically. This may lead to a loss of convergence or to 
instability, and it is preferable to evolve <f> with 



P 

-<- 

r 



where 



p = 2l + n. 



(8) 



(9) 



Therefore p is an even integer in an odd number (in par- 
ticular, three) of space dimensions, and an odd integer in 
an even number of space dimensions. With the first-order 
auxiliary variables 



i> = 4> 



this is equivalent to 



tj) = 7T , TT = ifj' + p— . 

r 



Regularity implies 

n(-r, t) = vr(r, t), ip(-r, t) = -t/)(r, t), 



(10) 



(11) 



(12) 



if we formally extend the functions to negative values of 
r. 

In this paper, we use the summation by part ideas 
(SBP) of Strand Q to construct finite differencing 
schemes of a given accuracy for (jl 1 [) . for positive integer 
p, with the regularity conditions (|12l) at r = 0, which ad- 
mit a discrete energy whose time derivative is a boundary 
term at a spherical outer boundary r = R. We then use 
the projection method of Olsson [2( to impose boundary 
conditions of one of the three forms 

pn + aip = 0, r = R, pcr>0, (13) 

(maximally dissipative boundary conditions), 

pn + [m' = 0, r = R, pp > 0, (14) 



or 



v (y + £v 



o. 



R, au > 0. 



(15) 



Any of these guarantee that a discrete energy exists that 
is positive definite and non-increasing. Hence the numer- 
ical scheme is stable. (A continuum energy exists and 
implies well-posedness also for the more general class of 
boundary conditions p7r + aip + [ip' + u(tp' + pip/r) = 0, 
but we have not been able to find an equivalent discrete 
energy.) 

We are aware of two other SBP schemes for (JTTJ) in the 
literature. 

The scheme proposed in Q and used in Q is SBP in- 
cluding at the outer boundary r — R, and is second-order 
accurate at fixed r. However, it has error terms that be- 
have as h 2 /r 2 . Therefore, at a fixed grid point i near 
the origin the error scales as 1/i 2 , and is independent of 
h. Therefore, the error does not converge in the maxi- 
mum norm, although for p ^ 2, 3 it converges with h 2 in 



the energy norm (|21|) defined below [5[. By contrast, the 
schemes we present here have an error of 0(h 2 ) or 0(h A ) 
at all grid points including the origin. 

A second-order accurate (including at the origin) finite 
differencing scheme for ([TTT) was described by Evans @ 
and is widely used, see for example [7[ for p = 2 and 
[H for p > 2. Although this method was not explicitly 
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designed as an SBP scheme, it actually is SBP at interior 
points and at the origin. This may explain its empirical 
robustness. Here, we complete it as an SBP scheme on 
the finite range < r < R by a suitable modification at 
the outer boundary. 

Our paper is structured as follows. Sec. |TT] discusses 
general considerations. In Sec. Ill Al we introduce our no- 
tation for the discretisation in r. In Sec. Ill Bl we identify 
the summation by parts (SBP) property that our differ- 
ence operators must obey for a discrete energy to exist. 
In Sec. Ill Cl we discuss the accuracy of these operators and 
in Sees . Ill Dl and UTEl we discuss the symmetry boundary 
r = and physical outer boundary r = R. 

In Sec. Mil we construct specific SBP schemes. In 
Sec. IIII Al we review the Evans method and make it SBP 
at the outer boundary. The Evans method is simple, but 
we have not been able to generalise it directly to 4-th or- 
der accuracy. Therefore, based on the general principles 
of Sec. m in Sec. MI Bl we derive another 2nd-order ac- 
curate SBP method, and in Sec. MI Cl we generalise this 
to 4th-order accuracy. It should be possible to construct 
implementations with arbitrary order of accuracy at in- 
terior points by extending our two examples. Sec. MI Dl 
confirms the expected accuracy and discrete energy con- 
servation of the three methods in numerical experiments. 

Sec. IIVI summarises our results and comments on their 
significance. 

Appendix [A] contains a rigorous discussion of ghost 
points at the origin, which are handled more informally 
in the main text. For completeness, Appendix iBl reviews 
Olsson's projection method [||. Appendix [C] uses the en- 
ergy method to show that a class of homogeneous bound- 
ary conditions that include first space derivatives (or first 
time derivatives, modulo the evolution equations), is sta- 
ble. In Appendix iDl we use the projection method to con- 
struct a discretisation for the three subclasses of bound- 
ary conditions stated above. 



II. GENERAL METHODS 
A. Discretisation 

Throughout this discussion we finite difference in r 
only, but assume the continuum limit in time. A fully dis- 
crete scheme is obtained by using a suitable ODE solver 
in t (the method of lines) . 

We use grid functions ^i(t) and Ui(t) on a grid rj to 
represent the continuum functions 7t(r, t) and ijj(r,t), as- 
suming that Ilj(i) = n(ri,t) and &i(t) — ip(ri,t), and 
that n(r, t) and ip(r,t) admit Taylor expansions in r to 
the required order at any r. From now on, we use a ma- 
trix notation where grid functions are written as column 
vectors, and finite differencing operators as matrices. We 
also suppress the t-dependence as it is relevant only later 
when we add time discretisation using the method of 
lines, that is we write w(r) and IL,, etc. 



A (2K + l)-point difference operator D is defined by 

i+s+K 

(M)< = £ DyVj (16) 

j—i+s — K 

where — K < s < K is an offset. The parameters 

of the difference operator are simply the elements of the 

band-diagonal matrix D. 

We assume a uniform grid with step size Ar = h. Our 
methods will require a grid that is either staggered about 
r = or centred about r = 0. In either case we find it 
convenient to introduce the notation 

1 3 

n=ih, i = M or * = 0,1,...M (17) 

that is, the grid index i takes half-integer values for the 
staggered grid and integer values for the centred grid. In 
either case tm = R- 

B. Summation by parts 

The continuum equations admit the energy 




with time derivative 

d § = [r p ^]t (19) 

where the second equation is obtained after using the 
evolution equations and integrating by parts. In order 
to control the growth of E, the boundary term at r = 
R must be controlled by a suitable boundary condition. 
The boundary condition 

pix -\- ci/j — 0, r = R (20) 

clearly guarantees dE/dt < for p, a > 0. Such bound- 
ary conditions are called maximally dissipative. We show 
in AppendixOthat; the more general boundary condition 
(IClj) also guarantees dEf, jdt < for a generalised energy 
E b . 

By contrast, the boundary at r = is not a physical 
boundary, but rather a point in the interior of the ball 
r < R. Indeed regularity alone implies that ip(0,t) = 0, 
and so the boundary term at r = vanishes. 

As is well-known, the continuum equations are well- 
posed in the norm provided by E because E is conserved. 
A summation by parts (SBP) finite differencing scheme 
exactly conserves a discrete equivalent E of the contin- 
uum energy E. This guarantees that it is stable (the 
discrete equivalent of well-posed) . 

We consider the discrete energy 

E = ^h p+1 (n^n + (21) 
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where 

W l = W, W l = W, W > 0, W > 0, (22) 

and write the finite differencing scheme as 

* = h^DU, II = h~ l D^. (23) 

The powers of h have been introduced so that W, W, D, 
D are all dimensionless and independent of h. The quan- 
tity r.i/h — i also has this property. The SBP property 
that guarantees that E is constant up to boundary terms 
is 

WD + {WD) 1 = B, (24) 
where the boundary operator B is defined by 



suggests the finite differencing operator 0] 



(i + ljPgj+i - (t - l) p ^-i 
2iP 



(31) 



However, for ^ = r and p = 2, the local error of the finite 
differencing operator h~ x D is exactly h 2 /r 2 — 1/i 2 . (For 
higher p > 0, terms up to (h/r) p also appear.) This does 
not go to zero with h at fixed i. 

We assume in the following that the operator h~ x D is 
a standard representation of d/dr with known accuracy, 
and concentrate on the difference operator D, where the 
p/r term may cause accuracy problems at the origin. 

The point r = 0, which arises (only) on a centered grid, 
must be treated specially. Taking the limit r —> of ((29)) 
at finite h we see that at r = the accuracy conditions 
are 



11*5* = X M p n A/ * 



M - 



(25) 



and the constant x obeys % — > 1 in the continuum limit 
M — > oo as h — > at fixed r = R. As is positive defi- 
nite, it is invertible, and we can consider D as determined 
by D, W, W and B: 



D = -W- x D l W + W^B. 



C. Accuracy 



(26) 



A useful reparameterisation of the difference operator 
coefficients D t j is obtained by expressing tyj — ip(rj) 
through the Taylor expansion of ip(r) about the fixed 
grid point r^. We can then write 

h- 1 (DV)i = coih-^in) + cwip'in) + c 2 ihip"(n) 

c 2 K,i^ 2K \n) + o(h 2K p7) 



■2K-1 



with a straighforward linear relation between the c a i and 
the Dij at each point i. In the following we adopt a 
simplified notation where the c a i (with a — 0, . . . ,2K) 
are written as c Q , ip(ri) simply as ip, etc., and fj simply 
as r. That is, we do not write the dependence on i, and 
all continuum quantities are evaluated at r — r{. 

Our difference operators are defined to be accurate to 
order 27V if they obey (using our abbreviated notation) 



2JV\ 



br x {DR)i = tt' + o(h 
h- 1 {b^) l = ^j + ^ 1 + o(h 2N ) 



(28) 
(29) 



A key observation for our problem is that the O (/i 2JV ) 
error needs to be obtained formally in the limit h — > 0, 
both at (approximately) constant r, and at constant i. 
The possible problem with the latter limit are error terms 
of the form h m /r n , which are 0(h m ) at constant r, but 
only 0(h m ~ n ) at constant i. As an example, the identity 



4' + -V 



(30) 



Ci = 1 +p, C 3 = C 5 = ... = C27V-1 = 0, 



(32) 



with the even c a undetermined because ip = "<!>" = • • • = 
for an odd function at the origin. For the methods with 
a centred grid that we have considered, these conditions 
will be discussed below for each specific method. 

Turning now to generic points with r > 0, naively one 
would expect the accuracy requirement for D to be equiv- 
alent to the following constraints on the coefficients of the 
difference operator (as defined above): 



ph 

co = — , ci = 1, c 2 
r 



C 2 N = 0. 



(33) 



Clearly, we need a stencil of width 2N + 1 or larger to 
control all these c a . 

However, we shall now see that we can replace some of 
the equalities (|3"3"|) with inequalities, which gives us ad- 
ditional freedom to impose the SBP conditions. Rather 
than devising a general notation, we present the cases 
N = 1 and N = 2. 

For N = 1, we make the following ansatz: 



ph fh 

co = h Oo - 

r V r 



ci = 1 - o" 



C2 = 01 



(34) 
(35) 
(36) 



where the 8 a may depend on i. The special case Sq = 
5i = gives (|3"3"|) . Substituting this ansatz into (|27|) 
gives 



r 2 (t-ip> 
r 



Here 



-8 x h 2 [r-V"] +i? 2 . 



R 2 = c 3 h 2 ^'" + c 4 h 3 ^"" + ..., 



(37) 
(38) 
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where for a 3-point stencil C3, C4, ...are known linear 
functions of cq, c\ and C2. Now, because of (fTHj) and 
regularity, both square brackets in (|37| are actually O(l) 
as r — > 0. Therefore, as long as So and S\ are bounded 
uniformly in i, the coefficients of h 2 in (|37[) are bounded. 
Similarly, as 03,04,... are regular functions of So and 5\, 
h 2 and all higher powers of h in ((38)) are also explicitly 
regular at r = and so we have the desired second-order 
accuracy. 

For N = 2 we make the ansatz 



Co 



Cl 



r V r 



I -S 



C3 



C4 



* 1 * 

r 



(39) 
(40) 
(41) 
(42) 
(43) 



which gives 



+S h 4 



+<5i/i 4 



r 2 r 3 



+a 2 /i* ir-V"! +i? 4 



(44) 



where R4 = 0(h 4 ) in the sense discussed above. Again, 
all the square brackets are regular at r = 0, and so we 
have fourth order accuracy if and only if the 8 a i are 
bounded uniformly in i. 

The values of Si in the methods we construct below are 
plotted against % and p in Figs. [1] and [5] 

It is clear that this method extends to arbitrary N, 
giving N equations to be solved and N + 1 inequalities 
(uniform in i bounds on the S a i) to be then verified. 

Informally, our method can be described as "trading 
r for h" . It works because the terms in square brackets 
above are all 0(1) as r — > 0, which in turn requires ip(r, t) 
to be a regular odd function of r. 

Note that for r ^> h, our c a given by (|34M36I) and 
(I39H43[) approach the values (|33|) . Similarly, R 2 and R4 
approach the values they take for the standard stencils 
D given below in ((62)) and (f93|) . which are precisely the 
minimum width stencils for a given accuracy that obey 

133). 



D. 



The symmetry boundary r = 



In numerical simulations using polar coordinates one 
is faced with the fact that r = is a boundary of the 



numerical grid, but is not in fact a boundary of the phys- 
ical domain. As a result, there are (typically) no physi- 
cal boundary conditions one can or must impose in the 
continuum limit, but the numerical simulation does re- 
quire boundary conditions. These are derived from the 
assumption that the desired solution is not less diffcrcn- 
tiable at r = than for r > (usually implicitly assumed 
smooth). The standard general approach to such "sym- 
metry boundary conditions" or "regularity conditions" is 
to identify variables which under the assumption of reg- 
ularity are either even functions of r (for example 7r) or 
odd functions of r (for example ip), and then to extend 
standard centered finite difference methods to the bound- 
ary by extending the numerical grid into a small number 
of "ghost points" representing negative r and which are 
populated by the assumed even or odd parity of the grid 
functions. 

From a strict SBP point of view, there are no ghost 
points, and finite difference operators are necessarily 
skewed near the boundary. The fact that r = is not 
a physical boundary is represented by the fact that B is 
zero at the boundary r = 0. 

However, we find that the use of ghost points allows a 
clearer derivation of our results, in that we do not need 
to discuss r = explicitly as a boundary. Rather than 
introduce a few ghostpoints, for our derivation we ex- 
tend all grid objects from 1/2 or 0, ... , M to —M, . . . M, 
corresponding to — R < r < R. 

We extend the grid functions to negative i as 



IF; 



(45) 



Because of how W and W define E, we can assume with- 
out loss of generality that 



W-i.-j 



0, 



i,3 > 0, 

(46) 

and similarly for W. B is extended by B_m,-m = 
—Bmm- In Appendix IA1 we prove from these assump- 
tions that ((45)) holds at all times if and only if 



(47) 



When coding our method, we implement D and D with 
a few ghost points. Equivalently, the ghost points can be 
explicitly eliminated. A rigorous treatment is relegated 
to Appendix El as it introduces additional notation not 
required for our main argument. 



E. The outer boundary r = R 

We begin by reviewing the one-dimensional wave equa- 
tion 



7T = 1p , 1p — 7T , a < x < b 



(48) 



with energy 



(n 2 +yj 2 )dx, -^ = W>t (49) 
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[This corresponds to the case p = of (fTTj) ]. It is natural 
to discretize this symmetrically in tt and ip, that is 



with energy 



n = D *, * = d il 



s = - (n*w n + **w *) 



and SBP condition 



W a D + (W Do) 4 = B Q 



(50) 



(51) 



(52) 



Notice that this problem is translation-invariant, and so 
.Do and Wo will naturally be translation-invariant in the 
interior, except for finite-sized end blocks. Do and Wo 
with various orders of accuracy in the interior and at the 
boundaries have been constructed by Strand 

In (fTTj) with p > additional problems result because 
the equations are not translation-invariant but depend 
explicitly on r. In the previous two sections we have 
addressed these problems at interior points and at the 
pseudo-boundary r = 0. Intuitively, the problem of im- 
posing SBP at the outer boundary r — R should be simi- 
lar to those solved by Strand because with increasing res- 
olution r varies less over a finite number of near-boundary 
grid points, and we should be able to make use of Strand's 
results. We now formalise this idea. 

Strand provides, among others, norms Wo that are the 
unit matrix except near the boundaries. We shall later 
construct norms W and W' for the problem without an 
outer boundary (on the range < r < oo) which are di- 
agonal everywhere except near r — 0. In this subsection 
we denote their restrictions to < r < R by the same 
names. We now set 



W = W W', 
w = W W', 
B = B W', 
D = D , 



(53) 
(54) 
(55) 
(56) 



where Wo, Do and Bo are the Strand operators with di- 
agonal W on -R < r < R. Then from ((2BJ and ^ we 
obtain 



b = w'^DqW 1 . 



(57) 



Let us assume we have already shown that this is accurate 
to some given order at the interior points where Do takes 
a simple translation-invariant form. Let us also assume 
that Do is accurate to some order r at the boundary 
points, and let us assume that 



W' = diag \i p (1 + 0(r r ))] 



(58) 



near the boundary and similarly for W' . Then near the 
boundary, at r < R, for fixed R, 

D = Do + ^z [1 + 0(i- T )] = D + | [1 + 0(h- T )] , (59) 

which is accurate to order r. 



III. SPECIFIC EXAMPLES 
A. The "Evans method" 

1. Interior points 

We concentrate at first on interior grid points by con- 
sidering the range < r < oo, Furthermore, points near 
the boundary r = can be treated like interior points if 
we use ghost points, as discussed in Sec. Ill Dl We use the 
notation D' etc. for operators on — oo < r < oo, and D 
etc. for operators on < r < R. 

The identity 



— + - I V 

dr r 



(p + 1) 



d{rP^) 
d{rP+ 1 ) 



(60) 



motivates the difference operator 



h- 1 (£>'¥)< = (p + 1) 



-P+1 



(61) 



which is the generalisation to arbitrary even p of a dif- 
ference operator introduced in Q. It is second-order ac- 
curate including at r = 0. We combine it with the usual 
second-order accurate 3-point symmetric difference oper- 
ator 



IL; 



(DTI), = 



and IL for negative i are defined by 
now 



(62) 



Consider 



W' = diag(tcj), W' = dmg{vi). (63) 
The SBP formula (gSJ then gives 



(D'tf)i = 



2wi 



(64) 



Comparing with (|61D we see that the Evans method 
is of this form with (using our convention = ih) 



(i + 1)p +1 -(i- 1)p +1 
2(p+l) ' 



(65) 
(66) 



We note that these are well defined for all i including 
i = 0. Indeed, the accuracy conditions at the origin 
(|3"2"|) for a method with diagonal energy (|53^) reduce to 
v% = (1 +p)wo, which is easily seen to hold for the ansatz 
(|65I66|) for even p. Note that the Evans method does 
not work for odd p (the wave equation in even space 
dimensions) on a staggered grid, as then wo = 0. 
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2. The outer boundary 

We now address the outer boundary. Following the 
general prescription of Sec. Ill El we set 



/ 



D = 



\ 



2 " 2 

-- - 

2 u 2 



(67) 



1 1/ 



r 





■ 









V 


«M / 




diag 


. . .,W M -2, 


WM 


diag 


, • ■ ■ ,BM-2, 


VM- 



B 



W M 
VM \ 

2 J' 

Note that this gives x = 1. From (l24l) we read off 



(68) 

(69) 
(70) 



2WM-2 2WM-2 

VM-2 r, U M 



V 

We read off 



2w M - 



2w M -l 
vm-i vm 
wm 



(71) 



W M W M 

[l + 0(M -1 )] hip' 

+ [1 + 0(M- 1 )] ^ + 0(^ 2 ) 



= h 



V/ + 1^ + ° {h) 



(72) 



where the last equality holds because M = R/h, and 
hence this method is first-order accurate at the boundary 
point i = M. 



3. Generalization to higher accuracy? 

The identity ((60)) suggests a generalization of the Evans 
method to accuracy order 27V, discretizing d/dr +p/r as 



(p + 1) 



D{rP+ 1 )' 



(73) 



where D is some discretization of d/dr of accuracy order 
2N. If we take the norms 



Vi = l H , Wi 



1 



D(i p+1 ), 

p+1 V ; ' 



(74) 



then ([73")) also obeys the SBP property. However, for the 
minimal- width centered stencils D of order larger than 2 
this does not work. To see this we differentiate ip(r) = 



ar + br 3 using the operators D of accuracy N = 1 and 
N = 2. The discretization errors are, respectively, 



b(p + 3)(2p+ l)h 2 
2b{p + 2,)p(p-l){2p-l) h A 
15 r 2 



(75) 
(76) 



In the latter case we see an error of the form /i 4 /r 2 , which 
becomes h 2 near the centre. 

We have not been able to generalize Evan's method to 
avoid this type of singular error term. Now we present 
a new, general approach to this problem, which allows 
high-order SBP methods that are regular at the centre. 



B. Another 2nd-order accurate implementation 
(SBP2) 

1. General considerations 

There are many ways of solving both the SBP condi- 
tions and the accuracy conditions for any required order 
of accuracy. We choose as our priority to make D and 
D as narrow as possible. For given D, wider W and W 
make D [given by (|26p ] wider, so we try to make them as 
narrow as possible. Finally, we choose W to be diagonal 
in order to make it easy to invert. We note that an over- 
all constant factor in W (and its inverse in W) cancels 
out from D. 

Given that h~ l D represents d/dr, it is natural to make 
the extended finite difference operator D' the "standard" 
minimal width operator for a given order of accuracy 27V 
at interior points. 



2. Interior points 

We concentrate at first on interior grid points. We 
assume again that D' is given by (|52^). and we look for 
W and W' of the form 



where W—i 



W' = W' = diag(w l ), 

The SBP formula gives 



2wi 



(77) 



(78) 



where for negative j is defined by (|4"5|) . 

We can now read off cq , c\ and C2 in terms of Wi . The 
one equality contained in (|34ti36|) . namely 



(9 



-) Co + Cl = l 



■P, (79) 

gives the linear recurrence relation of degree 2 for Wi 

(i + l)w i+x - (t - l)wi-! = 2{p + l) Wi . (80) 

The other two accuracy conditions are taken to define 
Si and 82 in terms of «;,-. On a staggered grid, we have 
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w_i/2 — ifi/2- We initially fix an arbitary value for W1/2, 
and can then solve the recursion for Wi for all i > 3/2. 
(Note that D' is unchanged if W is multiplied a constant 
factor) . On a centered grid, evaluating the accuracy con- 
ditions (|3"2"j) with = wi gives u>i = (1 + p)wq. We 
initially fix an arbitrary value of wq and can then solve 
the recursion for Wi for all z > 2. We note in passing 
that the method [3f is given at interior points by w; = i p , 
which does not obey ([50)1 . 

To work with a bounded quantity, we define the new 
variable Wi as 



wt = i p Wi . 



(81) 



(Therefore W-i = —Wi for odd p, while «;,• > 0.) It obeys 
the linear recurrence relation 



W; 



2(p + l) 



1-- 



p+i 

I - - ) „>;_ 2 . (82) 



Trying asymptotic solutions of the form 

00 

= ^ E C ™i k ~ m (83) 



m=0 



for constants p and fc shows that the two linearly inde- 
pendent solutions have p = ±1 and are (fixing a constant 
overall factor, and assuming i > 0) 



4 (+) = i+^4^+^-), 



12z 



(84) 



, } (-) = ( _i)i r 2(p+i) h _ (p+l)(p + 2)(p + 3) 

+o(r*)V 



(85) 



The first one is asymptotically constant, and the second is 
an oscillating decaying solution. The general asymptotic 
solution is an arbitrary linear combination of those, and 
hence it is also asymptotically constant. [The asymptoti- 
cally constant mode (|84[) for given p is a finite polynomial 



in i 2 of (the integer part of) l+p/2 terms. For example, 
restricting to i > 0, for p = 1 we have — 1 and for 

p = 2 we have w\ +) = 1 + l/(2i 2 ).] 

From these results we can infer the asymptotic be- 
haviour of the 8 a i. If Wi tends to a constant 7^ 
then we have 



= 7, 1" °{l )) 



I + 0(z- 2 ) 



(86) 
(87) 



Only in the case that u>oo = and only the oscillating 
mode is present is there a divergence in <5i, namely 



5oi = 2i 2 + 
x P + 2 



(p + 2)(p + 3) 
2 

+ 0^). 



0{i~ 



(88) 
(89) 



However, with our initial data tu_i/2 = w 1 j 2 on the stag- 
gered grid or wi = (1 + p)wq on the centred grid, the 
constant solution is present, and hence the sequence S a i 
converges as i — > 00, and is therefore bounded. Further- 
more, the upper bound of its absolute value is close to 
the asymptotic value, as we show in Fig. [T] 

Finally, we adjust the arbitrary overall factor such that 
lim.i_i.oo Wi = 1. On the centred grid we need 



p\ 



(90) 



for any value of p. Wi for i > is then actually given by 
the asymptotically constant polynomial (For even 

p, this is true for all i, but not for i = with odd p, 
where the special form of the accuracy condition at the 
centre needs to be used.) 

On the staggered grid we need for even p 



Wl/2 



[(P+1)!!] 2 
p + 1 



(91) 



which also leads to the polynomials (1541 . However for 
odd p the symmetry condition at the centre is incompat- 
ible with having only the asymptotically constant mode, 
and we need a contribution from the oscillating mode 
(|85p. For lim^oo Wi — 1 we now need 



h/2 



2 [(P+1)H] 2 

7T p+1 



(92) 



3. The outer boundary 



This works exactly as in Sec. IIII A 21 except that Vi = 
Wi and x — vm — t»M - > 1 only in the continuum limit. 



C. A 4th-order accurate implementation (SBP4) 

1. Interior points 

For N = 2, we have found a scheme where D 1 is the 
unique 4-th order accurate 5-point difference operator, 

(DTD, = 8(n ' +1 ' n - l) 12 (n - +2 ' II - 2) , (93) 

and where W is diagonal and W' is band-diagonal with 
three bands. We parameterize them as 

Wi,i = Wi, w-i = Wi, (94) 

W l . i = v^i = Vi, (95) 

Wi,i + i = u l+1 / 2 , u-i = Ui, (96) 

= Ui-1/2, (97) 

and all other components zero, where on the staggered 
grid the index on v and w takes half-integer values and 



12 5 10 20 50 100 

i 



1 2 5 10 20 50 100 

i 




FIG. 1. Values of So and Si for our second-order accurate (N — 1) methods, for p = 1, . . . , 10. SBP2 is in the left column and 
Evans in the right column. The staggered grids (half-integer i) and centred grids (integer i) are shown on the same plot. In all 
cases increasing values of p correspond to increasing \Si\, with even values of p shown in blue (dark) and odd values in orange 
(light). Note that the Evans method does not exist for odd p on the centred grid, and the corresponding dots are absent. We 
see a rapid convergence towards the respective asymptotic values (|86[l and (|87[l for SBP2, and So — > (p + 2)p(l — p)/3/(p + 1) 
and Si — > p/2 for the Evans method. Note that the i axis is logarithmic. 



the index on it takes integer values, and the other way 
around on the centred grid. % = on the staggered 
grid and iti/ 2 = U-1/2 = on the centred grid because 
we assumed in ([4"rJ)) that the off-diagonal blocks of W 
vanish. 

If we take temporarily as given, we can solve the 
accuracy conditions for the S a i and Wi, plus a linear re- 
currence relation of order 4 for Uj. On the staggered 
grid, we can fix Vi/% = v_x/2 and u 3/2 = u -3/2 arbitrar- 
ily, and solve the recurrence relation for Vi for i > 5/2 
starting from those four points and our choice of Uj . On 
the centred grid, the accuracy conditions (|32|) at the ori- 
gin reduce to (1 + p)w = V\ - (l/8)u 3 / 2 + (5/8)u 5 / 2 
and v 2 = vi + (63/8)u 3 / 2 — (27/8)u 5 / 2 - We can fix V\ 
and choose the the Uj arbitrarily and then compute Vi 
for i > 3 from the recurrence relation. (Note that Vq 
and U1/2 multiply which vanishes, and hence do not 
participate in the recurrence.) 



As wc shall see, it is sufficient to set all 



except 



for u\ on the staggered grid and it 3 / 2 and m 5 / 2 on the 



centered grid. In compact form, we can write 

{&*)< = (98) 
12wi 

where we have introduced the shorthand 

*l/2 =«l/2*l/2+«l*3/2, (99) 

*3/2 =f 3 /2*3/2+Ul*l/2, (100) 

^i = Vi^i, z>5/2. (101) 
for the staggered grid and 

*o = 0, (102) 

*1 = +U 3 /2*2, (103) 

$ 2 si; 2 $ 2 + !i3 /2 $i+ti 5/2 *3, (104) 

* 3 = V 3 V 3 + U 5/2 ^2, (105) 

^i=Vi%, z>4. (106) 
for the centred grid. 
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With the equivalent of (J5TJ) and ([531) for v i} the 
fourth-order linear recurrence for Vi has four independent 
asymptotic solutions with 



(2p - l)(p - 5) (p - 3)(p - 2) (p - l)p(p + 1) 



P=h -1, 



15, 4 



(107) 



all with = 0. The linearity of the recurrence implies 
that the general solution Vi is a linear combination of 
the four corresponding modes vf . It is possible to show 
that if the linear combination contains any contribution 
of the growing or oscillating modes then the 5 a are not 
bounded. Hence we must find a solution which only 
contains the asymptotically constant and the decaying 
modes. The freedom in ui, V1/2 and v 3 / 2 on the stag- 
gered grid, and in u 3 / 2 , U5/2 and V\ on the centred grid 
allows us precisely to cancel simultaneously the growing 
mode and the oscillating non-decaying mode and fix an 
overall constant factor. To do that we proceed as follows. 

We first compute three arbitrary solutions of the recur- 
rence up to some high value of i, say 1000. For example, 
on the staggered grid we can set each of ui, V\/2 and d 3 / 2 
to 1 and the other two to 0. The three solutions are dom- 
inated by the growing mode, and reach very high values, 
of order (4+V15) 1000 ~ 10 896 . We have detected extreme 
sensitivy of the solution in the initial conditions, roughly 
losing one decimal digit of precision per iteration, and 
hence the recurrence is solved with exact rational arith- 
metics, using Mathematica. 

Then we compute the asymptotic form of the modes, 
up to order 0(i~ s ). For instance for the asymptotically 
constant mode we have 



,-,(1) 



1 



(2p-l)(p-l)p(p+l)(p + 3) 
60i 4 

(2p - 3) (p - 3) (p - 2) (p - l)p(p + 1) (p + 3) 



504i 6 



+ o(r 8 ). 



(108) 



(for p — 2 this is simply 1 + 3/ (2i 4 ) + 0(i~ s ).) However, 
in contrast to the N = 1 case, these are only finite poly- 
nomials for odd p, but not for even p. For i ~ 1000 this 
expression will give results correct up to relative errors 
smaller than 10~ 18 for p < 10. We take three such values 
of i and construct a linear system to find which linear 
combination of our three solutions gives that mode t> t - . 
For such high values of i we can neglect the contribu- 
tion of the decaying mode. In this way we determine the 
values of Vi up to i = 1000. For larger i, and p < 10, 
the asymptotic series are accurate to 16 digits. In our 
experiments below we shall use up to p = 22, for wich 
values up to i = 2000 must be computed to use the given 
asymptotic expansions with relative errors below double 
precision. Note that we do not know if these series are 
convergent. 

From Vi we can compute This gives the following 
asymptotic behaviour for Wi, 



1 



(2p+l)(p+l)p(p-l)(p-3) 
60i 4 



504i 6 



+ 0{i- s ). 



(109) 



For i > 9/2 on the staggered grid and i > 6 on the 
centred grid the 5 can be computed from m and Wi as 
follows, 

x .5 Vj-2 ~ 8^-1 + 8^+1 - V i+2 .4 , n . 

Ooi = 1 — P l 1 (HO) 



Su = i 



hi = t ■ 



I2wi 

2 Vi- 2 - Vj-1 — Vj+l + Vi+2 

9wi 

2vi^ 2 - Vi-i + v i+ i - 2v i+2 



36wi 



(111) 
(112) 



The previous expansions imply that the S are bounded 
and have finite limits: 



-p(p-lf + 0(i- 2 ) 



5u = P -^ + 0(r% 

s 2i = -| + o(i- 2 ). 

6 



(113) 
(114) 
(115) 



The limit value of 5q is cubic in p. That means that 5q is 
very large for large values of p. Comparing with N = 1 it 
is plausible that So has an asymptotic limit which grows 
like p N+1 . 

We provide in our webpage Q data files with double- 
precision results for v, w and the S for 1 << 22 and i < 
2000. Formulas (|108p and (|109|l can be used to compute 
v and w for these p and i > 2000 to 16 digits. 

We find that for p = 1 and p = 2 (the wave equation 
in cylindrical and spherical symmetry), on the staggered 
grid, V1/2 < 0, so that W is then not positive definite. 
This problem is absent for p > 3, or on the centred grid. 
It is possible that allowing for ui other than u\ to be 
nonzero, this could be fixed, but we have not tried this. 



2. The outer boundary 

Strand provides boundary operators with interior accu- 
racy order 2N (2s in the notation of Strand) and different 
boundary accuracy orders r, depending on the number r 
of boundary points used and the form of the Wq (W in 
the notation of Strand) operator. He finds that r < N 
with r ranging from r = r + 1 for a general Wq to r = 2t 
for diagonal Wo . In our case 27V = 4 and we want to have 
diagonal Wo, so there are two possibilities: (r = 1, r = 2) 
and (r = 2,r = 4). 

For t = 1, following the general prescription of 
Sec. |HEJ we set 



D = 



j_ 

12 



_2 o 2 _J_ 

_ jl 8 i 2 

13 ^3 _ U y 1 § 

5 5 5 



(116) 
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FIG. 2. Values of 5o, 5i and 8% for SBP4 on the centred grid (left column) and on the staggered grid (right column), with 
p — 1, ... ,10. In all cases increasing values of p correspond to lines further from the axis Si = 0, with even values of p shown 
in blue (dark) and odd values in orange (light). Again, we see a rapid convergence towards their respective asymptotic values 
(|113Hll5)l . The fact that u\ (for the staggered grid grid) and U3/2 or tt 5 / 2 (for the centred grid) appear explicitly in the 
recurrence for a few low i points, but not beyond, produces some irregular behaviour at those points. 



B = 



V 

W — diag 



(117) 



Vm I 

,w M - 2 ,^-,—j, (118) 



W = diag . . .,v M -2, 



l3v M -i 5^m 
12 ' 12 



(119) 
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From we read off 



D 



2f M-3 
3w M -2 
"11-3 

13tu A f-i 





&VM-2 

13w M -i 

VM-2 

5uim 



2v M -\ 

3w M ~2 



7'UAf-l 

5ia>ai 



7«m 



This D is fourth-order accurate in the interior, and first- 
order accurate at i = M — 1, M. 



13w M _i 

(120) 



For t = 2, following the general prescription we set 



W = diag(. . . , v M - 



A9v M -3 43w M -2 59v M -i 17«m- 



48 



48 



48 



48 



(121) 



2v M - 



D 



12wm 








3m>a 

4 I'M- 



49tUM-3 












32u M -4 

49«>Af _ 3 
4t)Af-4 
43-UJA-J-2 






2r, 



3t«A f -4 



59VM-3 

86w M -2 



3f A-/-3 

34«>Af 



12m>j 
59?; a^ 



98w M -3 



UAf-2 



2lUA 
■iv M 





59uAf-i 
86wm-2 



59v m-1 
3Awm 





3vm 
98w M -3 

4vm. 
43w M -2 

VM 

2-wm - 1 
24n A f 
17 w M 



(122) 



The expressions for D and W are obtained by setting 
Vi and Wi to 1. 

D. Numerical tests 

We have tested our 2nd and 4th order accurate imple- 
mentations as described above, combined with 4th-order 
Runge-Kutta discretisation in time, with a Courant fac- 
tor At/h = 1/4. In SBP4 we use the numerical coeffi- 
cients Dij, or equivalently u\, Vi and u>i, calculated by 
the relaxation method described above up to i ~ 2000, 
and using the asymptotic results (|108l) . (|109[) for larger 
i. 

Initial data are set with tt a smooth function of com- 
pact support bounded away from r = and R, and 
if) = Q, so that at early times we can verify accuracy 
without the complication of boundary conditions. 

To determine the accuracy of the method we have in- 
vestigated self-convergence. This means that we check 
that e4 as defined by 

Ui(t) = v{r t ,t) + h A e A { ni t) + oih*) (123) 

(the Richardson expansion, here assuming 4th-order ac- 
curacy) is approximately independent of h. Starting from 
a reference resolution, we refine twice (by a factor of 2 on 
the centred grid, and a factor of 3 on the staggered grid). 
We can then obtain two independent estimates of e^ir, t) 
(or e3 or e2 if we are testing for third or second-order 
accuracy) and we compare these. Using refinement by 
a factor of 3 allows us to combine a staggered grid with 
fixed tm — R, and still have appropriate grid points of 



all refined grids align with the coarsest grid. Keeping 
R exactly resolution-independent is necessary to com- 
pare different resolutions, while aligned grids avoid the 
need for high-order interpolation to the same values r 
in comparing different resolutions. Both tt and ip and 
their associated errors scale as r~~ p l 2 to leading order, 
and we therefore plot r p / 2 Ti and r p / 2 e4. This is sensible 
also because technically speaking we expect convergence 
to the continuum in the energy norm E, rather than the 
usual L 2 norm, or in other words we expect convergence 
of r p / 2 7r and r p / 2 ip in the standard L 2 norm. 

We have tested Evans, SBP2, SBP41 and SBP42, for 
< p < 22, on the staggered and centred grids, with rep- 
resentative outer boundary conditions (see below). We 
note two exceptions. 1. The Evans method does not ex- 
ist on the centred grid for odd p, because wo is then not 
defined. 2. The energy of SBP4 is not positive definite 
on the staggered grid for p = 1,2, and so we would not 
expect it to be stable. However, we do not see signs of 
instability in our numerical experiments. 

In all other cases, our numerics confirm theoretical ex- 
pectations, as follows. At short times, before the solution 
has interacted with the boundary r = R, we see conver- 
gence at precisely the expected order, that is second or- 
der for Evans and SBP2, and fourth order for SBP41 and 
SBP42. As expected, the error is identical for SBP41 
and SBP42, because they are the same method except at 
the boundary. Surprisingly, the error e2 is almost identi- 
cal also between Evans and SBP2. This may be because 
they agree more and more with each other, and the naive 
method, with increasing i. 

Once the solution has interacted with the boundary, we 
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still see precise second-order convergence for Evans and 
SBP2. For SBP41 we see second-order convergence and 
for SBP42 we see third-order convergence. Note that in 
each case this global accuracy is one order higher than the 
accuracy at the boundary. We do not know if this was to 
be expected, see Sec. 12.7 of IjJ. For SBP41 and SBP42 
after the wave has interacted with the boundary, the error 
is dominated by grid noise, and hence the strict Richard- 
son expansion is lost, but the envelope of 63 or e2 respec- 
tively is still approximately resolution- independent. 

We have also evaluated the discrete energy at every 
time step. With the boundary conditions pix + fnr' = 
and aip + v(jp' + pip/r) = discussed in Appendix [Dl we 
have dEb/dt = 0. With the maximally dissipative bound- 
ary condition pit+aip = we have dE/dt — x^m^m < 0, 
and we have evolved the expected value of E by dis- 
cretising this in t using 4th-order Runge-Kutta. In all 
these cases the discrepancy between the evaluated and 
predicted numerical energy is of relative size 10 -8 , essen- 
tially independent of the method (Evans, SBP2, SBP41 
or SBP42) and of h, and increases linearly with t. As 
expected, it is due only to accumulated roundoff error 
(machine precision). 

IV. CONCLUSIONS 

The wave equations ([7]), © or (fTTj) appear in many 
areas of mathematical physics. In a larger class of ap- 
plications, these evolution equations are modified by r- 
dependent coefficients (for example in the wave equation 
on a curved spacetime), lower-order nonlinear terms (for 
example a scalar field potential), or couplings through 
lower-order terms to other equations in a larger system 
(for example stellar perturbation theory). 

Our experience suggests that the key stability prob- 
lems in such systems of equations arise already in the 
simple linear system considered here. Prior to the results 
presented here, we are not aware of any completely stable 
finite differencing scheme for these equations. Methods 
which appear to be stable at first fail to converge at suf- 
ficiently high resolution or large p, with problems either 
at the origin r = or the outer boundary r = R. 

It is surprising that the lower order term pip/r in (fTTj) 
alone can make standard schemes unstable, and that an 
elaborate SBP scheme is necessary. Note however that 
a standard centered finite difference implementation of 
the one-dimensional wave equation is already SBP except 
possibly at the boundaries, while the equivalent naive 
finite differencing of (|lip for p > is not SBP even 
at interior points. It seems highly unlikely to us that 
any boundary conditions can be devised to make such a 
scheme stable. 

Finally, the reader may be used to adding numerical 
dissipation to keep numerical methods stable. However, 
strictly linear equations such as the ones we consider here 
should not require explicit dissipation, which should be 
reserved for dealing with nonlinearities. Moreover, in a 



naive finite differencing of our equations, the amount of 
dissipation required increases rapidly with p, so that this 
is not a feasible approach if large p are required. 

The answer is then two- fold: first, to make the numer- 
ical scheme SBP at interior points and second, to make it 
SBP also at the boundary and to use the Olsson projec- 
tion method to impose suitable boundary conditions. In 
hindsight, the Evans method is already SBP at interior 
points, which explains why it works much better than 
a naive scheme. However, SBP at the outer boundary 
is also necessary to make it into a demonstrably stable 
scheme. This is not difficult to achieve, but without an 
explicit reference to discrete energy methods would not 
be at all obvious, especially as SBP methods appear to 
require a reduction in the order accuracy at boundary 
points [1]. Furthermore, the Evans method does not have 
an obvious generalisation to higher order accuracy, and 
does not work for odd p on a centred grid. 

As the main result of this paper, we have therefore 
constructed 2nd and 4th-order accurate SBP methods 
by using discrete energies that are non-diagonal near the 
origin. These can be made 1st and 2nd-order accurate, 
respectively, at the boundary r = R, leading to global 
2nd and 3rd-order accuracy after the solution has inter- 
acted with the boundary. (Our SBP4 is not positive def- 
inite for p = 1, 2 on a staggered grid.) It is clear that our 
ansatz can be generalized to order 2N in the interior and 
order TV at the boundary, giving rise to global accuracy 
of order N + 1. 

In the Appendix, we have used the SBP property to 
construct discretisations of the boundary conditions ([TBI 
IT51) that give rise to a strictly non-increasing discrete en- 
ergy, and hence guarantee numerical stability. 
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Appendix A: Rigorous treatment of ghost points 

Consider for physical grid points i > 0. We can 
write the use of ghost points explicitly as 

= h- 1 £ (DijUj + A,-,n-,) = h- 1 Dpi,, 

j>0 j>0 

(Al) 

where 

IK, Hj ■ />, h *,j>0. (A2) 
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We think of this as "folding over the ghost points". A 
similar observation holds for D, except that as is odd, 
the equivalent of (IA2[) is 



(A3) 



Note that D^ ^ even for i,j > 0, thus requiring a 
separate symbol. (The symbol D + is a reminder of the 
range i,j > 0.) The split of Df- into D^ and -Di,-j for 
i,j>0 is in general not unique. We do, however, have a 
natural prescription for this split if we assume that 
is translation-invariant, i.e. depends only on i — j even 
at the boundary. 

In order to extend Djj and Dij to negative i, we use 
the requirement that (|45|) hold at all times, or 



n_i = it, 



(A4) 



The first equation of (|45|) and the second equation of (|A4|) 
immediately give, for i,j > 0, that 



l>,, I I),. s I I) ... • l> ,. . 



= 0. 



(A5) 



The second equation of (|4"5)) and the first equation of 
(IA4|) . after substituting ([2"B| and using (j4"o]l . give, for 
i,j > 0, that 



(PF" 1 )^ (D kl - D fe> _i - £_ fei + D_ fc ,_,) = 0, 



and hence 



(A6) 



(A7) 



Taking the sum and difference of (|A5|) and (|A7[) . we ob- 
tain (g7J). 

Finally, we have, for i,j > 0, that 

= J2 -(W'^ikiDik - D_ lk )Wij 

k,l>0 

= -(W~ 1 )ifc(Afc + A,-fc)W' w 

M>0 



^ -(W-^ifcp 4 -*)*!^ 



(A8) 



fej>0 



and so P4]l holds for the operators D + and D + with 
ghost points folded in if and only if it holds for the ex- 
tened operators D and D. This confirms that the intro- 
duction of ghost points is just a matter of notation. 



Appendix B: The projection method for imposing 
boundary conditions 

For completeness, this Appendix summarises relevant 
methods from 0. Suppose a first-order in space and 
time system of PDEs in one spatial dimension has been 
discretised in space as 



Note that the vector u in general ranges over multiple 
variables (for example it and ip) as well as grid points 
(for example i), and we use calligraphic letters such as V 
for operators on this vector space. 

Suppose this system has a discrete energy 



E = -utyVu 



and obeys the SBP property that 

B = i (WD + 2?*W) 

is a boundary operator. Then 

dE t 
—— = u Bu 
dt 



(B2) 



(B3) 



(B4) 



is a boundary term. 

We want to impose one or several homogenous linear 
boundary conditions that we write as 



Lu = 0. 



(B5) 



In matrix notation where a is a column vector, C is a 
matrix that has one row for each boundary condition. 

To focus on the essential idea, we define the inner prod- 
uct 

(u,v)=u t Wv. (B6) 

In this notation we can write 

1 dE 
E=-(u,u), — = (u,Vu). (B7) 

The adjoint with respect to this inner product is defined 
by 

(Au,v) = (u,A f v), (B8) 
and is therefore given by 

= W- l A l W. (B9) 

The operator 

V = l-W- 1 £ t (£W- 1 £ t )- 1 £ (B10) 
clearly obeys 

V 2 =V, CV = 0, V^=V, (Bll) 

and so is a self-adjoint projection operator into the space 
of grid functions that obey the boundary conditions. If 
we now use the semi-discrete evolution equation 



u = VDu 



(B12) 



u = Du 



(Bl) 



instead of (|B1[) , we have £ii = exactly, and hence £u — 
and therefore Vu = u at all times if it holds initially. 
Then we have 

dE 

— = (u, VVu) = (Vu, Vu) = (u, Vu) (B13) 

as before and so both the discrete energy bound and the 
desired boundary conditions hold. 



15 



Appendix C: Continuum boundary conditions 
involving derivatives 



Appendix D: Numerical boundary conditions 
involving derivatives 



Consider the class of boundary conditions of the form 
pir + atp + fnr' + v (V + ^) =0, r = R (CI) 



or equivalcntly 

piT + (Tip + pip + Z/7T 



0, 



R 



(C2) 



for p, (j, p, v not all vanishing at once. To fix an overall 
sign, we also assume that at least one of them is posi- 
tive. We now use an energy argument to show that these 
boundary conditions give rise to a stable initial-boundary 
value problem if p,a, p,v > with pp + av > 0. [The 
maximally dissipative special case p = v = with pa > 
is also stable based on the energy E defined in (|18|)]. 
We consider the energy 



R p 

E b = E+—(piP + vn) 2 r=R , 



(C3) 



where E is given by (|T%|). E b stands for E modified by a 
boundary term, and s is 



pp 



av. 



(C4) 



Its time derivative is 



dEb 
dt 



= R p 



R p 



1 



nip H — (pip + vk) ( Wljj + VTT 
s \ 

i 

vir) (pir + aip) 



r=R 



nip (pip - 

s 



r=R 



RP f ,<2 2 \ 

(palp + Vp-K ) r - 

S 



R- 



(C5) 



The necessary and sufficient conditions for Ef, to be pos- 
itive definite and its time derivative to be non-positive 
are 



p>0, cr>0, p>0, v>Q. 



(C6) 



We have dE b /dt — if pa = pir = 0, and dE b /dt < 
otherwise. However, the limiting case pp = av = is 
not allowed because it would give s = 0, except for the 
maximally dissipative sub-case p = v = 0, where E and 
dE/dt are given by (HSJ) and (O instead of (|U3]) and 
(IC51). 



We define the modified numerical energy 



Eh — E 



X h p M p 
2s 



(At* 



vll 



M 



(Dl) 



where x parameterises finite differencing error in the 
boundary term, as defined by ([25]) . We find 



dEb 
dt 



xhPMP 



n M *Af 



+i (/x* M + jTIm) (m^a/ + ^n M ) , (D2) 



so if the numerical boundary could be chosen to be 

pn M + cnpM + pipM + vkm = 0, (D3) 

the argument could be completed as in the continuum 
case. However, ii — VT>u and not T>u. We have not been 
able to find an ansatz for C and E b such that dEb/dt < 0. 

Consider however the two subclasses of boundary con- 
ditions where dEb/dt = in the continuum. Consider 
first the case a = v = with pp > 0. Then (|D2j) reduces 
to 



at \ p 



= X h p MP [ U m ^m + ^Mh-^D^M 
P 



(D4) 



where the second equality holds because in this special 
case the boundary condition is independent of ^ so that 
V only acts on n, and the last equality holds if we im- 
plement Cu = as 



pli 



A I 



ph-\DU)M = 0. 



(D5) 



The case p = p = with va > works the same way, 
with the roles of n and ^ interchanged. 
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